1. Field of the Invention
This invention relates generally to the field of geophysical prospecting. More particularly, the invention relates to the field of seismic data processing. Specifically, the invention is a method for imaging prestack seismic data for seismic migration.
2. Description of the Related Art
Seismic migration is the process of constructing the reflector surfaces defining the subterranean earth layers of interest from the recorded seismic data. Thus, the process of migration provides an image of the earth in depth or time. It is intended to account for both positioning (kinematic) and amplitude (dynamic) effects associated with the transmission and reflection of seismic energy from seismic sources to seismic receivers. Although vertical variations in velocity are the most common, lateral variations are also encountered. These velocity variations arise from several causes, including differential compaction of the earth layers, uplift on the flanks of salt diapers, and variation in depositional dynamics between proximal (shaly) and distal (sandy to carbonaceous) shelf locations.
It has been recognized for many years in the geophysical processing industry that seismic migration should be performed pre-stack and in the depth domain, rather than post-stack, in order to obtain optimal, stacked images in areas with structural complexity. In addition, pre-stack depth migration offers the advantage of optimally preparing the data for subsequent AVO (Amplitude Versus Offset) analysis. Pre-stack depth migration has traditionally been performed through the application of Kirchhoff methods. However, because of recent advances in computing hardware and improvements in the design of efficient extrapolators, methods that are based on solutions of the one-way wave equation are starting to be used.
It has been well established in the literature that wave equation-based methods offer the kinematics advantage of implicitly including multipathing, and can thus more accurately delineate structures underlying a complex overburden. However, there has been considerably less discussion of the dynamic advantages that might (or should) be realized through wave equation-based methods. This is not surprising, since Kirchhoff pre-stack depth migration has undergone a much longer evolution than its wave equation-based counterparts. Part of this evolution has been the development of various factors or strategies to account for geometrical divergence and spatial irregularities in acquisition.
Poststack migration is often inadequate for areas that are geologically complex, because the migration is performed at a late stage, after the seismic data have been affected by the NMO (normal moveout) or DMO (dip moveout) corrections and CDP (common depth point) stacking. One prestack option is depth migration of common shot gathers, or shot records. This method is thus called “shot record” or “shot profile” migration. Shot record migration can give more accurate imagining, better preserve dip, and provide accurate prestack amplitude information. These properties have made shot record migration one of the more popular methods of wave theory migration.
Shot record migration is described in Reshef, Moshe and Kosloff, Dan, “Migration of common shot gathers”, Geophysics, Vol. 51, No. 2 (February, 1986), pp. 324-331. Reshef and Kosloff (1986) describe three methods for migrating common-shot records. Discussion of the accuracy and efficiency of shot record migration compared to full pre-stack migration is contained in the two articles: Jean-Paul Jeannot, 1988, “Full prestack versus shot record migration: practical aspects”, SEG Expanded Abstracts, pp. 966-968 and C. P. A. Wapenaar and A. J. Berkhout, 1987, “Full prestack versus shot record migration”, SEG Expanded Abstracts.
Seismic migration comprises two steps: (1) wave extrapolation and (2) imaging. The present application deals with the second step of imaging in seismic migration. A corresponding method of wave extrapolation in seismic migration is described in co-pending U.S. patent application, “Method for Downward Extrapolation of Prestack Seismic Data”, by the same inventor and assigned to the same assignee as the present patent application.
The key elements necessary for an accurate shot record migration may be derived by examining the formal representation for the reflection response from a single, flat interface. Shot record migration is then seen to be the set of operations necessary to invert for the plane-wave reflectivity. After subsequent generalization to the multi-layered case, these elements include (1) definition of the source, (2) extrapolation of the wavefields' phases over depth with an amplitude-compensating algorithm that preserves vertical energy flux, and (3) application of an imaging principle that is based in some way on the definition of reflectivity.
Let the upper half-space contain a point, compressional seismic source with coordinates xs≡(xs,ys,zs) and spectrum S(ω), and let it also contain seismic receivers, such as hydrophones, with coordinates xr≡(xr,yr,zr). Depth is chosen to be positive downward, and for simplicity assume that zr=zs. Finally, let the depth of the interface be z′. A far-field expression for the pressure of the upgoing (reflected) wavefield, designated by U, is then given by:
                              U          ⁡                      (                                          x                r                            ,                              y                r                            ,                              z                s                            ,              t                        )                          ≈                              (                          i              /                              c                0                2                                      )                    ⁢                                    ∫                              -                ∞                            ∞                        ⁢                                                            S                  ⁡                                      (                    ω                    )                                                  ·                                  exp                  ⁡                                      (                                                                  -                        ⅈ                                            ⁢                                                                                          ⁢                      ω                      ⁢                                                                                          ⁢                      t                                        )                                                              ⁢                                                ∫                                      -                    ∞                                    ∞                                ⁢                                                      exp                    ⁡                                          [                                              ⅈ                        ⁢                                                                                                  ⁢                                                                              k                            x                                                    ⁡                                                      (                                                                                          x                                r                                                            -                                                              x                                s                                                                                      )                                                                                              ]                                                        ·                                                            ∫                                              -                        ∞                                            ∞                                        ⁢                                                                  exp                        ⁡                                                  [                                                      ⅈ                            ⁢                                                                                                                  ⁢                                                                                          k                                y                                                            ⁡                                                              (                                                                                                      y                                    r                                                                    -                                                                      y                                    s                                                                                                  )                                                                                                              ]                                                                    ⁢                                              R                        ⁡                                                  [                                                                                    exp                              ⁡                                                              (                                                                                                      -                                    2                                                                    ⁢                                  ⅈ                                  ⁢                                                                                                                                          ⁢                                  ω                                  ⁢                                                                                                                                          ⁢                                                                      ξ                                    0                                                                    ⁢                                                                                                                                                                                        z                                        ′                                                                            -                                                                              z                                        s                                                                                                                                                                                                                )                                                                                                                    ω                              ⁢                                                                                                                          ⁢                                                              ξ                                0                                                                                                              ]                                                                    ⁢                                              ⅆ                                                  k                          y                                                                    ⁢                                              ⅆ                                                  k                          x                                                                    ⁢                                              ⅆ                        ω                                                                                                                                                    (        1        )            where c0 and ξ0≡ξ(c0) are the velocity and vertical slowness in the upper half-space, ω is the radial frequency, kx and ky are horizontal wavenumbers, and R is the plane-wave reflectivity at the reflector. In an acoustic medium, reflectivity R has the simple form:
                              R          ⁡                      (                                                            k                  x                                ω                            ,                                                k                  y                                ω                            ,                              z                ′                                      )                          =                              (                                                            ρ                  1                                                  ξ                  1                                            -                                                ρ                  0                                                  ξ                  0                                                      )                    /                      (                                                            ρ                  1                                                  ξ                  1                                            +                                                ρ                  0                                                  ξ                  0                                                      )                                              (        2        )            where ρ0 is the density in the upper half-space and ρ1, c1, and ξ1≡ξ(c1) are the density, velocity, and vertical slowness, respectively, in the lower half-space. Likewise, the downgoing wavefield, designated by D, at the reflector depth z′ has the form:
                              D          ⁡                      (                                          k                x                            ,                              k                y                            ,              z              ,              ω                        )                          ≈                              (                          i              /                              c                0                2                                      )                    ⁢                      S            ⁡                          (              ω              )                                ⁢                      exp            ⁡                          [                              -                                  ⅈ                  ⁡                                      (                                                                                            k                          x                                                ⁢                                                  x                          s                                                                    +                                                                        k                          y                                                ⁢                                                  y                          s                                                                                      )                                                              ]                                ⁢                                    exp              ⁡                              (                                                      -                    ⅈ                                    ⁢                                                                          ⁢                  ω                  ⁢                                                                          ⁢                                      ξ                    0                                    ⁢                  z                                )                                                    ω              ⁢                                                          ⁢                              ξ                0                                                                        (        3        )            
Equation (3) provides a representation that is convenient for theoretical analysis, but it is not the most practical form for numerical application. Since extrapolation for shot record migration is commonly performed in the (ω,x) domain, it is more convenient to represent the source in the same domain. Away from the source, the wavefield is given by the inverse transformation of Equation (3):
                              D          ⁡                      (                                          x                ′                            ,              ω                        )                          =                              S            ⁡                          (              ω              )                                ⁢                                    exp              ⁡                              (                                  ⅈ                  ⁢                                                                          ⁢                  ω                  ⁢                                                                                                                                    x                          ′                                                -                                                  x                          s                                                                                                            /                                          c                      0                                                                      )                                                                    c                0                2                            ⁢                                                                                    x                    ′                                    -                                      x                    s                                                                                                                          (        4        )            where x′=(x′,y′,z′) specifies the observation point for the wavefield in the subsurface.
Equation (4) can be used in at least two different ways. First, the wavefield can be specified at the bottom of a homogeneous overburden and then downward continued through the underlying, heterogeneous structure. Alternatively, in a land setting, a wavefield can be specified at some shallow depth using the local velocity, reverse-extrapolated back to the surface and then downward continued through the underlying, heterogeneous medium.
Thus, the steps needed to invert the recorded seismic data for the plane-wave reflectivity, R, can be identified. First, the source wavefield is modeled to the reflector depth, perhaps through Equation (4). Second, the upgoing, recorded wavefield (U) is extrapolated to the reflector depth z′. Third, the reflectivity R is computed by
                              R          ⁡                      (                                          k                x                            ,                              k                y                            ,                              z                ′                            ,              ω                        )                          =                              U            ⁡                          (                                                k                  x                                ,                                  k                  y                                ,                                  z                  ′                                ,                ω                            )                                            D            ⁡                          (                                                k                  x                                ,                                  k                  y                                ,                                  z                  ′                                ,                ω                            )                                                          (        5        )            and then inverse transformed back to the spatial domain.
Several imaging principles for shot record migration have been reported in the literature. Five common methods are reviewed. Kinematic imaging principles only require the upgoing, receiver wavefield U, while dynamic imaging principles depend on some functional relationship between both the upgoing and downgoing, source wavefields U and D.
A popular, kinematic approach is given by:
                                          R            ⁡                          (              x              )                                =                                    1              n                        ⁢                                          ∑                                  j                  =                  1                                n                            ⁢                                                U                  ⁡                                      (                                          x                      ,                                              ω                        j                                                              )                                                  ⁢                                  exp                  ⁡                                      [                                          ⅈ                      ⁢                                                                                          ⁢                                              ω                        j                                            ⁢                                              τ                        ⁡                                                  (                          x                          )                                                                                      ]                                                                                      ,                            (        6        )            where x=(x,y,z) correponds to the image point and τ(x) is the one-way time from the source to the image point as determined from a ray trace. One can expect it to yield results that are partially subject to the limitations of ray theory. In order to account for spreading effects, factors must be computed from the ray trace and applied in the same manner as for Kirchhoff migration. For a further description of this approach, see, for example, the two articles Reshef, Moshe and Kosloff, Dan, “Migration of common shot gathers”, Geophysics, Vol. 51, No. 2 (February, 1986), pp. 324-331, or Reshef, Moshe, “Prestack depth imaging of three-dimensional shot gathers”, Geophysics, Vol. 56, No. 8, (August, 1991), pp. 1158-1163.
An average correlation in the (x,ω) domain is given by:
                                          R            ⁡                          (              x              )                                =                                    1              n                        ⁢                                          ∑                                  j                  =                  1                                n                            ⁢                              [                                                      U                    ⁡                                          (                                              x                        ,                                                  ω                          j                                                                    )                                                        ⁢                                                            D                      *                                        ⁡                                          (                                              x                        ,                                                  ω                          j                                                                    )                                                                      ]                                                    ,                            (        7        )            where the superscript “*” denotes complex conjugation. This approach will yield poorly-resolved (low frequency), but robust images. However, it does not yield a valid measure of reflectivity. An accurate recovery of divergence effects requires that the reflectivity be computed as some form of a ratio between the source and receiver wavefields.
The classical definition of reflectivity in the (x,ω) domain is given by a deconvolution:
                              R          ⁡                      (            x            )                          =                              1            n                    ⁢                                    ∑                              j                =                1                            n                        ⁢                                                            U                  ⁡                                      (                                          x                      ,                                              ω                        j                                                              )                                                                    D                  ⁡                                      (                                          x                      ,                                              ω                        j                                                              )                                                              .                                                          (        8        )            This imaging condition is very close to the theoretically correct method, which entails a computation in the wavenumber, rather than the spatial, domain. This method can yield “floating point errors” where the source illumination is very weak.
A normalized deconvolution approach in the (x,ω) domain is given by:
                                          R            ⁡                          (              x              )                                =                                    1              n                        ⁢                                          ∑                                  j                  =                  1                                n                            ⁢                                                                    U                    ⁡                                          (                                              x                        ,                                                  ω                          j                                                                    )                                                        ⁢                                                            D                      *                                        ⁡                                          (                                              x                        ,                                                  ω                          j                                                                    )                                                                                                                              D                    ⁡                                          (                                              x                        ,                                                  ω                          j                                                                    )                                                                                                                            ,                            (        9        )            where |D|=(DD*)1/2. This algorithm is less sensitive to noise than the deconvolution approach in Equation (8).
After introducing pre-whitening, Equation (9) becomes:
                              R          ⁡                      (            x            )                          =                              1            n                    ⁢                                    ∑                              j                =                1                            n                        ⁢                                                                                U                    ⁡                                          (                                              x                        ,                                                  ω                          j                                                                    )                                                        ⁢                                                            D                      *                                        ⁡                                          (                                              x                        ,                                                  w                          j                                                                    )                                                                                                                                                                                D                        ⁡                                                  (                                                      x                            ,                                                          ω                              j                                                                                )                                                                                                            2                                    +                  ɛ                                            .                                                          (        10        )            Equation (10) is probably the most commonly used imaging method in the industry. The pre-whitening term ε stabilizes the solution in the presence of noise. However, it has the undesirable side effect of inhibiting an accurate recovery of geometrical divergence. The optimal value for ε is highly dependent on the data and varies with each shot record and also with spatial position. Considerable testing is thus required to find a single, “best compromise” value. It should be chosen as small as possible while still yielding a stable, clean-looking solution. For an overview of the approaches discussed with reference to Equations (7)-(10), see, for example, Jacobs, Allan, 1982, “The prestack migration of profiles”, Ph.D. thesis, Stanford Univ., SEP-34, and the classic book, Claerbout, Jon, “Imaging the Earth's Interior”, Blackwell Scientific Publications, Ltd., Oxford, United Kingdom, 1985, 1999.
Valenciano, Alejandro, Biondi, Biondo, and Guitton, Antoine, 2002, “Multidimensional imaging condition for shot profile migration”, Stanford Exploration Project, Report SEP-111, Jun. 10, 2002, pp. 71-81, formulate a type of least-squares inverse for the imaging condition in shot record migration. A more robust (smoothed) least-squares inverse might be formulated by considering data at spatially neighboring locations in computing the inverse at location x. However, their inverses are all performed in the time domain, rather than the frequency domain. In the frequency domain, convolution and deconvolution correspond to multiplication and division, respectively. In the time domain, deconvolution involves the construction of matrices with time-shifted operators. The least-squares inverse formed in Valenciano et al. (2002) is based on these matrices.
Valenciano et al. (2002) also propose a simple variation of the pre-whitening factor used in Equation (10) above. A binary (0. or 1.) weight can be added to eliminate ε when the product U·D* exceeds some threshold and retain it when U·D* is unacceptably small. Results in areas of clear imaging would thus not be biased by the pre-whitening, and “floating point errors” would be avoided at other locations.
Thus, a need exists for a method for imaging prestack seismic data in shot record migration that can be fully automated and recover amplitudes in heterogeneous media.